Two-dimensional algorithm of the density matrix renormalization group 
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We propose a new approach to implement the density matrix renormalization group (DMRG) in 
two dimensions. With this approach the initial blocks of a L x L lattice are built up directly from the 
matrix elements of a (L — 1) X {L — 1) lattice and the topological characteristics of two dimensional 
lattices is preserved in the iteration of DMRG. By applying it to the spin-1/2 Heisenberg model 
on both square and triangle lattices, we find that this approach is significantly more efficient and 
accurate than other two-dimensional DMRG methods currently in use. 



I. INTRODUCTION 

The density matrix renormalization group (DMRG) is 
an optimized iterative numerical method. Since its devel- 
opment by White in 1992 pL this method has achieved 
tremendous success in studying ground state properties 
of one-dimensional (ID) interacting electrons. It has also 
been successfully extended to finite temperatures , to 
momentum space 0, and to the calculation of dynamic 
correlation functions 

The DMRG starts from a small system which can be 
handled rigourously. A large chain, called superblock, is 
then built up from this small system by adding a num- 
ber of sites at a time. At each stage, the superblock 
consists of system and environment blocks in addition to 
a number of extra sites. Graphically, a superblock can be 
represented as {S » s * e E), where S and E represent the 
system and environment blocks and » s and » e the extra 
sites added to S and E, respectively. S and » s (similarly 
E and » e ) form an augmented block, which becomes the 
system (environment) block in the next iteration. How- 
ever, in order to keep the size of the superblock basis from 
growing, the basis for the augmented blocks is truncated. 
Hence the DMRG is a basis truncation method. However, 
unlike the conventional renormalization group method, 
the truncation is done for each augmented subblock and 
the basis states retained are determined not by their ener- 
gies but by their probabilities projecting onto the ground 
state (or other targeted states) of the superblock. These 
probabilities are determined by the reduced density ma- 
trix of the augmented system (or environment) block. 

To construct the density matrix, the ground state \ip) 
of the superblock is first diagonalized with the Lanczos 
or other sparse matrix diagonalization algorithm. The 
reduced density matrix of the augmented system (or en- 
vironment) is defined by tracing out from \tp) (ip\ all the 
degrees of freedom that do not belong to this block: 



P = Tr (i 
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(1) 



Thus (E © » e ) is considered as a statistical bath to the 
augmented system. The density matrix is semi-positive 
definite. Rs eigenvalue is equal to the projection proba- 
bility of the corresponding eigenvector in i.e. 



A l = X;i(Ai,e J #>| S 



(2) 



where (A;, |A;)) is an eigenpair of p and {|e 3 )} is a basis 
set of (E © » e ). 

Given the density matrix, an entropy can be defined 
for the augmented system according to the standard ther- 
modynamic relation 



S = —Trp In p = — A; In A; 



(3) 



The maximum of the function /(A) = —A In A is located 
at A = e _1 . When < A < e _1 , /(A) increases mono- 



tonically with A. When A > 



/(A) decreases with 



A. No more than two A; can be larger than e since 
A; = 1. Thus if the contribution to the entropy from 
the largest A; is larger than that from the largest dis- 
carded eigenvalue of p, the DMRG is also a maximum 
entropy method. 

There are two approaches in forming a superblock. In 
literature they are often referred as the finite and infinite 
lattice approaches. In the infinite lattice approach in one 
dimension, the environment block is generally chosen as 
the space reflection of the system. In the finite lattice 
approach, the size of the superblock is fixed and the en- 
vironment block is chosen as the remaining part of the 
lattice for a given system block. The infinite lattice ap- 
proach allows the size of the superblock to be flexible and 
can be used to study the thermodynamic limit directly. 
However, the finite size approach is more accurate in cal- 
culating quantities for a system with fixed lattice size. 

The DMRG can also be used to study thermodynamic 
properties of a ID quantum || or 2D classical system 
|||. In this case, the transfer matrix of a Hamiltonian 
system, instead of the Hamiltonian itself, is diagonalized. 
The free energy and other thermodynamic quantities are 
determined by the maximum eigenvalue of the transfer 
matrix. The transfer-matrix DMRG method treats di- 
rectly an infinite lattice system and has therefore no finite 
lattice size effect. 

A simple extension of the DMRG to more than ID 
would be to replace the single sites added between the 
blocks with a row of sites, either along a principal axis 
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H or along a diagonal [||. However, the extra degrees 
of freedom added to the system would make the size of 
the Hilbert space prohibitively large. Therefore, the two- 
dimensional algorithm should be developed so that only 
a single site is added to each subblock at a time. 

In practice the extension of the DMRG to more than 
ID is to map a higher dimensional lattice onto a ID one, 
namely to choose a path to order all lattice sites Jl(]] . The 
mapping breaks the lattice symmetry and introduces long 
range interactions among lattice sites. Therefore, the 2D 
procedure differs from the ID one in that there are addi- 
tional connections between the system and environment 
blocks. 

A typical mapping, as illustrated in Fig. 1, is to fold 
a ID zipper into 2D. This is basically a multi-chain ap- 
proach since the length of the folded zipper is unlimited 
but the width is fixed. For a 2D gas of non-interacting 
electrons, Liang and Pang found that the number of 
states needed to maintain a certain accuracy grows expo- 
nentially with the width of the lattice . This conver- 
gence was also confirmed for an algorithm where a row of 
sites was added at each step || . Although no proof has 
been given, this statement is often referred to as most 
probably valid for any 2D DMRG calculation. 



This leads to a strong restriction on the basis states and 
allows the number of states kept to increase substantially. 
Unlike its real space counterpart, the momentum space 
DMRG treats the kinetic energy rigorously. Hence this 
method works better in the weak coupling limit. How- 
ever, the application of the momentum space DMRG has 
its own limitations. For example, it is very difficult, if 
not completely impossible, to apply this method to a 
pure spin system like the Heisenberg model. 

In this paper, we introduce a new approach to imple- 
ment the DMRG in real space in 2D. Instead of ordering 
the lattice sites row by row as in the multi-chain ap- 
proach, we order the lattice sites by the order along the 
diagonal direction. As will be shown later, this is a truly 
two-dimensional method which allows us to build up the 
initial system and environment of a L x L lattice system 
based on the results on a (L — 1) x (L — 1) lattice and is 
particularly suitable for handling 2D lattice models. 

The rest of the paper is arranged as the following. In 
Sec. II a truly 2D algorithm of the DMRG is introduced. 
In Sec. Ill, as an example of the application of the algo- 
rithm, the ground state energy of the spin- 1/2 Heisenberg 
model is evaluated on both square and triangle lattices. 
The study is summarized in Sec. IV. 
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FIG. 1. A superblock in a "multichain" algorithm. The 
system and environment blocks are enclosed by dashed lines. 
Black spots are the sites added to the system and environ- 
ment. 

This multi-chain approach is simple to implement in 
the DMRG iteration. However, with this approach, the 
calculations on (L — 1) x (L — 1) and L x L are per- 
formed independently. The information obtained from 
the iterations on a (L — 1) x (L — 1) lattice is not used 
in the preparation of the initial sub-block matrices in the 
calculation for a L x L lattice. This is undoubtedly a 
loss of the efficiency. It may result in the loss of the 
accuracy as well, since the topological characteristics of 
square lattices is not well manifested in the preparation 
of the initial block states and the sweeping procedure of 
DMRG iterations. 

The momentum space DMRG provides an alternative 
way to implement the DMRG in two or higher dimensions 
W. In this representation the momentum is conserved. 



II. A 2D ALGORITHM OF THE DMRG 

In this section we will take the square lattice as an ex- 
ample to show how to build up initial blocks of a L x L 
lattice from a (L — 1) X (L — 1) lattice. The extension to 
any 2D lattice which can be topologically transformed to 
a square lattice by adding or removing some of the near- 
est or next nearest neighbor interactions from the square 
lattice, such as triangle, hexagonal and Kagomi lattices 
(an example for such a transformation is given in Fig. ||), 
is straightforward. 
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FIG. 2. A L x L (L = 5 here) triangle lattice (a) can be 
taken as a L x L square lattice with extra next-nearest neigh- 
bor coupling (b). 

Let us start from a 2 x 2 lattice. Fig. ^a shows the 
order of the sites after the 2D — > ID mapping. As the 
system is small, the Hamiltonian can be fully diagonal- 
ized. 
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FIG. 3. (a) is a 2 x 2 lattice, (b) is the initial configura- 
tion of superblock for a 3 x 3 lattice system. As indicated by 
the numbers, the lattice sites are ordered along the diagonal 
direction. The initial system contains three sites linked by 
the solid line in the lower left corner. The initial environment 
contains four sites, also linked by the solid line, at the upper 
right corner, (c) same as for (b) but for the next iteration. 
Black spots are the extra sites added into the superblock. 

Fig. ||(b) shows the configuration of the initial su- 
perblock for a 3 x 3 lattice system. As indicated by the 
number shown in the figure, the lattice sites are ordered 
from the lower left corner to the uppper right corner 
along the diagonal. The initial system contains three 
sites linked by the solid line in the lower left corner. The 
initial environment contains all the four sites in the up- 
per right 2x2 lattice. All the matrix elements for these 
initial subblocks can be obtained from the results previ- 
ously obtained on the 2x2 lattice. We add site 4 to the 
system and site 6 to the environment to form the aug- 
mented system and environment blocks. Unlike in a real 
ID system, these two added sites are not nearest neigh- 
bors in the mapped ID system. After a standard DMRG 
calculation for this superblock, the augmented system 
block can be updated and taken as the new system in 
the next iteration. 

In the next iteration (Fig. ||c) the system contains four 
sites (i.e. sites 1-4) and the environment contains only 
three sites at the upper right corner (i.e. sites 7-9). Since 
the two sites (i.e. sites 5 and 6) to be added to the system 
and environment are nearest neighbors in the mapped ID 
lattice, from now on the DMRG finite system sweeping 
can be done exactly as in a true ID system. 

Similarly, the DMRG iterations on a 4 x 4 lattice can 
be done based on the results of the 3x3 lattice. As for 
a 3 x 3 lattice, a 4 x 4 lattice (Fig. ||a) can be formed by 
two corner cut off 3 x 3 lattices with two isolated sites. 
The initial system contains 6 sites linked by a solid line 
in the lower left corner (i.e. sites 1 — 6) and the initial en- 
vironment contains 8 sites, also linked by a solid line, in 
the upper right 3x3 lattice (i.e. sites 8, 9, 11 — 16). The 
configurations of these two blocks can be found from the 
previously studied 3x3 lattice with or without a space 
reflection. We add site 7 to the system and site 10 to 
the environment to form the augmented system and en- 
vironment blocks. Again, these two sites are not nearest 



neighbors in the mapped ID system. But the standard 
DMRG calculation can be done as usual. The augmented 
system block is then updated and taken as the new sys- 
tem in the next iteration. 




(a) (b) 



FIG. 4. (a) a 4 x 4 lattice decomposed as two partially over- 
lapped 3x3 lattices (enclosed by the short dashed squares) 
and two sites at the two corners outsides these 3x3 lattices 
(i.e. sites 6 and 10). The number besides each lattice site 
gives the order in the mapped ID system. The sites in a 
system or environment block are linked by solid lines. Black 
spots are the sites added in. (b) same as for (a) but for the 
next iteration. The environment (sites 10-16) is a space re- 
flection of the system (sites 1-7) with respect to the center of 
the 4x4 lattice. 

In the next iteration (Fig. |^b), the augmented system 
in the last iteration becomes the new system. It contains 
seven sites (i.e. sites 1 — 7). In this case, since the total 
number of sites in the environment is also seven, the en- 
vironment can therefore be taken as the space reflection 
of the system with respect to the center of the 4x4 lat- 
tice, i.e. sites 10 — 16. All the matrix elements of this 
environment can be obtained from the space reflection of 
the system. The sites now added into the system and 
environment are the two nearest neighboring sites in the 
mapped ID system. Thus starting from this iteration, 
the standard finite system sweeping can be done as in a 
ID system, without considering how the 4x4 lattice is 
constructed from the 3x3 lattices. 

The above procedure can be repeated to larger square 
lattices. In general, the initial superblocks of a L x L 
lattice can be formed based on the results of the system 
and environment blocks in a (L — 1) x [L — 1) lattice. 
We order all the lattice sites like a folded zipper with 
inequal width along the diagonal. If the first site at the 
lower left corner of the L x L lattice is labeled as 1, then 
the two sites to be added in will have the coordinates 
Xi = (L-l)L/2 + l andX 2 = L(L + l)/2in the mapped 
ID system, respectively. (An example is given in Fig. ^ 
for a 5 x 5 lattice system.) We take the first (L — l)L/2 
sites in the lower left corner as the initial system and all 
the sites in the upper right (L — 1) x (L— 1) square lattice 
not used by the system as the initial environment. The 
DMRG calculation can be done as before. The system is 
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always augmented and updated. At the first few itera- 
tions, the site which is added to the environment is hxed 
at X2 and is not exactly next to X\ in the mapped ID 
lattice. This continues until the environment can be gen- 
erated by the center reflection of the system and the two 
sites added to these two blocks become nearest neighbors 
in the mapped ID system. After that the standard finite 
system sweeping can be done as in an ordinary ID lattice. 




FIG. 5. The first three initial configurations of superblocks 
for a 5 x 5 lattice system. At the first two iterations, (a) and 
(b), the two sites which are added to the system and envi- 
ronement (i.e. X\ and X2) are not next to each other in the 
mapped ID system. At the third iteration, (c), X 1 and Xi 
become next to each other in the mapped ID system. 

III. THE 2D HEISENBERG MODEL 

In this section, we take the spin- 1/2 Heisenberg model 
as an example to demonstrate how good our approach 
is compared with the multi-chain approach. The ground 
state energies on both square and triangle lattices are 



evaluated. For these 2D systems, there are currently 
rather precise results available, mainly from large-scale 
Monte Carlo calculations and series expansions. There- 
fore the accuracy of our results can be assessed by com- 
paring with these results. 

The Heisenberg model is defined by the Hamiltonian 

# = E S *' S , (4) 

where Si is the spin operator and the summation runs 
over all nearest neighbors. In real space at the same pa- 
rameters and number of states, the truncation error in 
a system with periodic boundary conditions is usually 
much higher than with free boundary conditions, there- 
fore we use free boundary conditions. 

The total spin S 2 is a good quantum number for the 
isotropic Heisenberg model. This symmetry has been 
used in obtaining all the results presented below. We 
have also performed finite system iterations using both 
our algorithm and the multichain one. In the multichain 
calculations, we have used an algorithm introduced in 
Ref ]Tl[ | to build up the initial system or environment 
blocks. 

Table I compares the ground state energy per bond 
obtained by the true 2D approach, Eid, with that ob- 
tained by the multichain approach, E mci on both square 
and triangle lattices. For square lattices, E 2 d is always 
lower than E mc . Since the DMRG satisfies the varia- 
tional principle, this means that the true 2D results are 
more accurate than the multichain ones. Moreover, the 
difference {E mc — Eid) / \E mc \ increases with increasing 
lattice. Thus the improvement of the true 2D approach 
over the multichain approach becomes more and more 
significant as the lattice size is increased. For triangle 
lattices, Eid is slightly higher than E mc when L is small. 
However, for large lattices Eid is much more accurate 
than E mc . The increase of (E mc — Eid) I \E m c\ with in- 
creasing size in the triangle lattice is even larger than in 
the square one. 

For a given L, an accurate estimate of the ground state 
energy (similarly other physical quantities) can be ob- 
tained by extrapolating Eid to the limit m — > 00. This 
can also be done (l2) by extrapolating Eid with respect 
to the truncation error Ae, since the limit m — > 00 is 
equivalent to the limit Ae — > 0. The extrapolation with 
respect to the number of retained states is difficult to 
implement since the asymptotic behavior of Eid in the 
limit m — > 00 is unknown and there is some uncertainty 
in determining the function used in the extrapolation. 
However, we find that the Ae dependence of Eid is gen- 
erally very simple and can be well described by a power 
law in the limit Ae — ► 0. An example is given in Fig. ^ 
where the Ae dependence of Eid on a 6 x 6 square lattice 
is shown. In the figure, the solid line is a polynomial fit 
(up to the quadratic term in Ae) to the data. From the 
fit the ground state energy per bond for this 6x6 system 
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is estimated to be —0.36212. For other cases, this fitting 
procedure can be similarly done. 
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FIG. 6. The ground state energy per bond E211 as a func- 
tion of the truncation error for the spin- 1/2 Heisenberg model 
on a 6 x 6 square lattice with free boundary conditions. The 
solid line is a polynomial fit to the data. 

To obtain the ground state energy in the thermody- 
namic limit, we need to do a finite size scaling for the 
results obtained from the above extrapolation. In a pe- 
riodic system, the leading size correction to the ground 
state energy per bond is of order 1/L 3 . |l3) However, in 
an open system as considered here, the finite size effect is 
stronger and the leading size correction is of order 1/L. 




FIG. 7. Ground state energy E2d versus 1/L of the Heisen- 
berg model with free boundary conditions on square lattices. 
The behavior of Eid on even lattices is different to that on 
odd lattices. But the extrapolated value in the limit 1/L — ► 
is the same within numerical errors. 

Figs. (^) and (||) show the scaling behavior of 
the ground state energy on square and triangle lat- 



tices, respectively. For the square lattice, the extrap- 
olated ground state energy in the limit 1/L — > is 
Eoo w —0.3346. This agrees very well with the prob- 
ably best currently available estimate, obtained from 
large-scale quantum Monte Carlo calculations, of Eoo w 
-0.334719(3). Q The result of spin wave theory is 
Eoo = —0.33475 up to the fourth order correction 
p"5| . For the triangle lattice, the extrapolated result is 
Eoo ~ —0.1814. It is also consistent with the quantum 
Monte Carlo results obtained by Capriotti et al Jl6| , 
Eoo ~ -0.1819, and by Bernu et al g7|, Eoo ~ -0.1825. 
The second order spin wave result is Eoo = —0.1822. flSl] 
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FIG. 8. -©2d versus 1/L for the Heisenberg model with free 
boundary conditions on triangle lattices. The solid line is a 
polynomial fit to the data with mod(L, 3) = 0. 

The above comparison indicates that accurate results 
for the ground state energy can be obtained using the 
algorithm introduced above. In obtaining these results, 
the symmetry of the total spin S 2 is considered and up to 
300 states are retained. This calculation can be readily 
done on a moderate workstation. With the aid of modern 
parallel computers, we should be able to keep more num- 
ber of states (e.g. 3000 states) and to further increase 
the accuracy. 



IV. CONCLUSION 

We have developed a new approach to implement the 
real space DMRG in 2D. We point out that a L x L 
lattice can be taken as an assembly of two partially over- 
lapped (L — 1) x (L — 1) lattices plus two extra sites 
and therefore the initial blocks of a L x L system can be 
built up directly from the blocks of a (L — 1) x (L — 1) 
system. This is a truly 2D algorithm of the DMRG. It 
preserves a higher degree of the symmetry of 2D lattice 
than the multichain approach and can be readily used 
in the DMRG calculation. For the spin-1/2 Heisenberg 
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model on both square and triangle lattices, the ground 
state energies obtained with this approach are consistent 
with the quantum Monte Carlo results and better than 
those obtained with the multichain approach for large 
lattice systems. 
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TABLE I. Comparison of the ground state energy per bond of the Heisenberg model on square and triangle lattices with 
free boundary conditions obtained by the true 2D approach, E-zd, with that obtained by the multichain approach, E mc . m is 
the number of states retained. The lattice size is L? . 
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